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A three-dimensional compressible CFD code has been developed to study the heat transfer characteristics 
of a twin-power piston y-type Stirling engine. The results include temperature contours, velocity vectors, 
and distributions of local heat flux along solid boundaries at several important time steps as well as var¬ 
iation of average temperatures, integrated rates of heat input, heat output, and engine power. It is found 
that Impingement is the major heat transfer mechanism in the expansion and compression chambers, 
and the temperature distribution is highly non-uniform across the engine volume at any given moment. 
This study sheds light into the complex heat transfer mechanism inside the Stirling engine and is very 
helpful to the understanding of the fundamental process of the engine cycle. 

© 2014 Elsevier Ltd. All rights reserved. 


1. Introduction 

Stirling engine, a nearly 200 years old invention, has attracted 
much attention and been the subject of intensive research by many 
industrial and academic institutes in recent decades due to its 
many specific advantages [1-4], These advantages are: very high 
efficiency, silence in operation, low level of toxic emission if pow¬ 
ered by fuel, and the ability to use almost any kind of heat sources, 
including solar, biological, geothermal, or even industrial waste 
heat. Especially the last advantage makes it regarded as an impor¬ 
tant green energy device and potentially one of the most promising 
solutions to the global warming problem. However, there are also 
disadvantages associated with the Stirling engine. It is low in 
power/weight ratio, unable to change torque and power swiftly, 
and expensive (especially for those associated with concentration 
solar energy); and these have preclude the engine’s widespread 
into some of today’s important applications such as transportation 
and electricity generation. Despite these disadvantages, Stirling en¬ 
gine technology has been constantly improved to produce new 
generation of engines with higher power/weight ratio and higher 
power output and efficiency; and it might become more popular 
and widespread into more applications in the future. 

Improvement on Stirling engine’s performance relies heavily on 
extensive numerical and experimental research on engine cycle. To 
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design a new Stirling engine, a numerical analysis often precedes 
experiment to find some general directions to follow. It is also a 
more suitable tool than experiment to carry out parameter optimi¬ 
zation. Therefore, many numerical models for Stirling engine cycle 
have been developed. A detailed overview of the existing numeri¬ 
cal models for the analysis on different Stirling engine cycles can 
be found in Mahkamov [5], These models are basically classified 
as the first-, second-, and third-order models, according to the 
model hierarchy. Here, the term "order” has nothing to do with 
the mathematical terms in the governing equations. In these mod¬ 
els, the entire engine is divided into 3, or 5, or more sections (con¬ 
trol volumes), so they are practically zero- or one-dimensional 
models. Among them, the model proposed by Schmidt is the ste¬ 
reotype, and improved variants of Schmidt’s model are nowadays 
the most widely adopted models [6-8] for numerical analysis of 
Stirling engines. Some second- or third-order models take non-iso- 
thermal and transient effects into account and are capable of pre¬ 
dicting the average values of volume, pressure, and heat transfer 
rate in each section at each time step of the cycle, allowing the 
overall heat input and power output in a cycle to be calculated. 
However, because these models are zero- or one-dimensional 
models, they, at best, only resolve spatial variations in the axial 
direction of the engine. In this case, any variation in the transverse 
direction that is important for predicting the engine’s performance 
accurately will not be accounted for. Neither can they resolve the 
effects caused by changes in some geometrical features such as 
the shapes and dimensions of displacer and displacer cylinder, 
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Nomenclature 



Cp 

constant pressure specific heat (kj kg -1 K -1 ) 

t 

time (s) 

k 

thermal conductivity (Witt 1 K _1 ) 


time period of an engine cycle 

l, 

length of the power piston linkage bar (m) 

T 

temperature (I<) 

h 

length of the power piston connection rod (m) 

U, V, w 

velocity components respectively in x-, y-, and z-direc- 

h 

length of the displacer linkage bar (m) 


tions (ms -1 ) 

U 

length of the displacer connection rod (m) 

ii , z>, w 

local frame moving velocity components respectively in 

Ic 

height of displacer (m) 


x-, y-, and z-directions (ms -1 ) 

m 

mass (kg) 

V 

volume (m 3 ) 

p 

pressure (Pa) 

w 

engine power (W) 

Q 

heat transfer rate (W) 

x, y, z 

components of Cartesian coordinate system (m) 

Q 

heat flux or heat transfer rate per unit area (Wnr 2 ) 



R 

gas constant (Jkg -1 IC 1 ) 

Greeks 


R, 

outer radius of the power piston (m) 

P 

crank angle of the displacer 

r 2 

inner radius of the displacer cylinder (m) 

P 

density (kg m 3 ) 

Rd 

outer radius of the displacer (m) 

n 

engine efficiency 

Ri 

thermal conductivity ratio between solid and gas mate- 

0 

crank angle 


rials 

P 

viscosity (Pas) 

r, 

crank radius of the power piston (m) 

CD 

engine speed (rads -1 ) 

r 2 

crank radius of the displacer (m) 




regeneration chamber, or the internal gas circuits, all of which can 
only be defined in multi dimensions. This means that they are not 
very accurate in terms of the optimization of the geometries and 
dimensions of individual components during the design of a Stir¬ 
ling engine. In addition, heat transfer rates inside the expansion 
and compression chambers are often estimated by using constant 
convective heat transfer coefficients, a practice that has over-sim¬ 
plified the complexity of the real heat transfer mechanism taking 
place inside these chambers. These factors together with other 
simplifications result in poor predictive accuracy returned by these 
models. They tend to over-predict engine’s output power and effi¬ 
ciency by the margins from 30% to 100%, or even more. 

A Stirling cycle involves very complicated heat and mass trans¬ 
fer processes, and the engine itself consists of many multi-dimen¬ 
sional components whose geometrical effects can be important. 
None of these can be resolved in great details or very accurately 
by those aforementioned numerical models. Yet, resolving these 
sophisticated physical processes and the effects from some geo¬ 
metrical features are crucial to the prediction of the overall engine 
performance. Especially understanding the effects of the geometri¬ 
cal features is of practical importance because to build a new Stir¬ 
ling engine, engineers have to determine the detail geometries and 
dimensions of all components in the engine. Therefore, more elab¬ 
orated models are needed to improve the accuracy of numerical 
analysis. Computational fluid dynamics (CFD) approach is one of 
such numerical models by far that are capable of simulating those 
complicated processes taking place in Stirling engine and hence 
returning more accurate prediction on the overall engine perfor¬ 
mance. It is also capable of resolving the effects caused by changes 
in small geometrical features and potentially be a great tool for the 
design and development process of a new Stirling engine. How¬ 
ever, there exist many challenges associated with applying a CFD 
approach on a Stirling cycle. Although a Stirling engine is a close 
system where precise boundaries can be defined, the problem is 
unsteady, multiphasic, and multi-scale in nature. Furthermore, 
the engine involves moving parts such as displacer and power pis¬ 
ton; hence the volume of computational domain is changing with 
time. These bring challenges in two aspects. In the aspect of coding, 
simulating the reciprocating motions of displacer and power piston 
requires the use of moving mesh techniques which introduce com¬ 
plex algorithm and raise the difficulty of coding. In the aspect flow 
physics, the gas flow is compressible; the heat transfer mechanism 


is conjugate heat transfer, involving solid and fluid media; and in 
most cases, the flow is turbulent; and these give rise to complex 
heat and mass transfer processes. To resolve these processes, it is 
necessary to use very fine grids and small time steps and solve 
many coupled governing equations at the same time, resulting in 
problems of convergence and long CPU hours. Due to these chal¬ 
lenges, many attempts to apply CFD approach on Stirling engines 
focused on separated components such as combustion chamber 
of the heater [9], There have been very limited reports on using 
CFD approach to study the full cycle of a Stirling engine. Mahka- 
mov [5] used both second-order approach and axisymmetric CFD 
approach to analyze the working process of a solar Stirling engine. 
They found significant differences in temperature data and engine 
performance obtained by both methods. The transient temperature 
variations by CFD approach are far from the harmonic distributions 
given by the second-order approach, and CFD model returns much 
more accurate output power than the second-order model. Mahka- 
mov 10] also reported a 3D CFD approach on a prototype biomass 
Stirling engine to improvement its performance. The engine was 
also analyzed by the first-order and second-order approaches. 
The CFD method identified some key factors that reduce the power 
of the engine: a geometrical feature that causes high level of 
hydraulic losses in the regenerator, and the entrapment of the 
gas in the pipe connecting two parts of the compression space 
and the dead volume introduced by the pipe. Such detailed analy¬ 
sis on the effects imposed by multi-dimensional geometrical fea¬ 
tures can only be resolved by a CFD approach. The study also 
highlighted the usefulness of adopting a CFD approach to identify 
and solve problems during the design of a Stirling engine. 

Among different types of Stirling engines, the y-type engine has 
many specific advantages as described in Kongtragool and Wongw- 
ises [11]. It is a low to medium temperature differential engine; 
therefore it can run on temperature difference from just tens of de¬ 
gree to a couple of hundreds of degree. This is a very important 
advantage because the engine can be powered by a wide range 
of energy sources, including low-grade sources such as solar, bio¬ 
fuel, geothermal, or even industrial waste heat. It is also much eas¬ 
ier and cheaper to be constructed and maintained than other types 
of Stirling engines. Therefore, the y-type engine is very suitable for 
domestic use, which can potentially become a very large market, 
because of the readily availability and low cost of low-grade fuel 
and the engine’s simplicity to construct and maintain. Although 
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CFD approach has been applied to the research and development of 
industrial Stirling engines, in open literature, there has not been a 
detailed CFD study on the basic heat transfer processes in the y- 
type Stirling engine cycle. Yet, understanding the fundamental pro¬ 
cesses of its engine cycle in great detail is a key step to improve the 
engine’s performance. The objective of this study is to perform CFD 
simulation on the working cycle of a y-type Stirling engine to re¬ 
veal the characteristics of its heat transfer processes in great detail. 
This is an important step to further improve the performance of 
such engine cycle. 


2. Mathematical model 


The Stirling engine analyzed in this study is the twin-power pis¬ 
ton y-type Stirling engine reported in Chen, et al. [12,13], where 
the engine was experimentally studied and numerically analyzed 
by an improved variant of the second-order model. The schematic 
of the geometrical configuration of the engine is given in Fig. 1. 
According to Fig. 1, the displacements of the piston and displacer 
are: 

z p (9 ) = l c i - (-u sin 9+\Jl\- r] cos 2 9 + l 2 Y (1) 


Zd(jS) = /cl - 




( 2 ) 


where fi = 0 - 90°. The velocities of piston and displacer are: 


w p (9) = rococos9 - 


r^co cos 0 sin 0 
cos 2 9 


(3) 



Fig. 1. The configuration and definition of geometric parameters of the y-type 
Stirling engine. 


Wd(fS) = r 2 co cos fi 


rococos/? sin/I 



(4) 


The configuration of the computational domain, expressed by 
the computational mesh, is illustrated in Fig. 2. As seen, due to geo¬ 
metrical symmetry, only a quarter of the real engine space is re¬ 
quired. Since the engine in Chen’s [13] experiment is quite large, 
the flow inside the engine could become turbulent at relatively 
low engine speed. Therefore in this study, the engine has been 
scaled down by a half to reduce Reynolds number and retain lam¬ 
inar flow inside the engine, precluding the need to involve turbu¬ 
lence modeling. In addition, due to its smaller size, the engine is 
not equipped with a regenerator, and regeneration is performed 
by the working gas flowing forwards and backwards between the 
expansion and compression chambers through the narrow gap 
formed by the displacer and the inner wall of displacer cylinder; 
that is, the solid material of the wall of displacer cylinder is acting 
like regenerator material to store and release heat. Furthermore, 
due to the fact that the heat capacity of solid wall is much larger 
than that of the gas, the temperature of solid wall is further as¬ 
sumed to retain at a constant value (on the hot- or cold-end wall) 
or a fixed profile along the z-direction (on the displacer cylinder 
wall), hence eliminating the need to include the solid domain. 
Some other assumptions to facilitate CFD simulation are made as 
follows: fluid viscosity and thermal properties all materials are as¬ 
sumed constant, the working gas is air and assumed to be ideal gas, 
and effects from mechanical friction and thermal radiation are ig¬ 
nored. Therefore, the flow and heat transfer phenomena of the en¬ 
gine cycle can be solved by transient three-dimensional laminar 
and compressible Navier-Stokes equations plus energy equation 
and ideal gas equation of state as follows: 

Continuity equation : 

dp d d d 

i + lTx {pU) + lK {pV) + lK {pW) = 0 - (5) 



Fig. 2. The computational mesh, showing the configuration of the computational 
domain. 
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Momentum equations: 


d(pu) 

dt 


d_ 

dx 


(piiu) + 


d_ 

dx 


(, pvu ) + 


d_ 

dx 


(pwu) 


dp (d 2 u cPu 9 2 u\ 

dx + P \W + dy 2 + d^)’ 


d(pv ) d . - d - d ~ . 
-DT + dx {puv) + dx {pvv) + ¥x {pwv) 
dp (d 2 V cfv d 2 v\ 

~ ~dy + P \W + dy 2+ ~d^)' 


( 6 ) 


( 7 ) 


d(pu) d . d - d _ 

T + ^ uw) + ^ w)+ ^ ww) 


dp 

dz 


-in-PS + V 


dx 2 dy 2 


cPw 

dz 2 


Energy equation: 

d(pT) d d . d . 

-ir + d~x ipuT)+ irx (pvT)+ irx (pwT) 


d 2 T d 2 T d 2 T 

W + dy 2+ ^ : 


Equation of state: 
P = pRT- 


( 8 ) 


( 9 ) 


( 10 ) 


In the above, u = u- u b , v = v - v b , w = w-w b are relative veloc¬ 
ity components between fluid and local moving frame that moves 
with tit, Vb, w b , respectively in x, y, and z directions. In Eq. (9), the 
viscous dissipation effect has not been included. As argued in Mah- 
kamov [10], this effect is insignificant compared with other heat 
transfer effects in Stirling engine, thus it can be neglected. The ini¬ 
tial conditions are: 


0 = 0°, u = v = w= 0, p = 101.0 kPa, T = T L . (11) 

Boundary conditions are: 

At hot end, z = 0: 

u=v = w = 0, T = T h . (12) 

At cold end, z = l cl - l c2 and on the wall of power cylinder: 

u=v = w = 0, T = T l . (13) 

The surface of displacer is assumed adiabatic, thus the condi¬ 
tions on displacer surfaces are: 

u=v = 0, vv — w d , g = 0. (14) 

where n is the direction normal to the wall of displacer. The surface 
of power piston is also assumed adiabatic, giving: 

u=v = 0, w — Wp, | = 0. (15) 

The temperature on the lateral wall of the displacer cylinder 
(the wall of regenerative channel) is assumed to maintain at a fixed 
linear profile as: 

u=v = w = 0, T = Ti + n Z , A T h -T l ). (16) 

(kl — k2) 


3. Numerical procedure 

The numerical procedure in this study is based on an in-house 
unstructured-mesh, fully collocated, finite-volume code, 
‘USTREAM’ developed by the corresponding author. This is the 
descendent of the structured-mesh, multi-block code of ‘STREAM’ 


[14]. In a Stirling engine, the actions of the displacer and the power 
pistons cause expansion and compression of the engine volume, 
and a moving mesh technique is needed to simulate such varia¬ 
tions in volume. There are two commonly used moving mesh tech¬ 
niques: one is adding or deleting cells to change the volume of the 
computational domain, and the other is to directly change the 
geometries of cells to expand or compress their volumes. In this 
study, the second strategy is adopted. The number of cells in the 
computational domain is fixed, and the expansion or compression 
of engine volume is performed by moving the local boundaries of 
cells to increase or decrease the volumes of the cells, and collec¬ 
tively, they can simulate the expansion or compression of the en¬ 
gine volume just like the action of an accordion. The moving of 
cell boundary creates local moving frame velocity components 
u,v. w. However, only the component w exists due to the fact that 
the motions of displacer and power piston are both along the z- 
direction. 

Another important issue is to simulate the compressible flow 
inside the engine. Since USTREAM is a pressure-based SIMPLEC 
[15 finite volume code, the pressure correction equation has to 
be modified to account for variable density in compressible flows. 
The pressure-velocity coupling scheme in Lebon et al. [16] is 
adopted here to perform this task. 


4. Results and discussion 

The dimensions of the simulated y-type Stirling engine are 
given in Table 1. The working gas is air with a gas constant of 


Table 1 

The values of the geometrical parameters of the Stirling engine. 


R, (m) 

0.0128 

R 2 (m) 

0.0400 

Rd (m) 

0.0390 

r, (m) 

0.0200 

r 2 (m) 

0.0125 

1, (m) 

0.0650 

h (m) 

0.0225 

h (m) 

0.0475 

U (m) 

0.0775 

Id (m) 

0.0740 

la (m) 

0.2125 

lez (m) 

0.1105 



Fig. 3. p ave -V diagram of an adiabatic cycle of the y-type Stirling engine for code 
validation purpose. 
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Fig. 4. Variations of total, expansion chamber, compression chamber volumes 
versus crank angle within an engine cycle. 


287.0 Jkg -1 1< ', and the initial crank angle and gas pressure are 0.0 
degree and 101.0 kPa, respectively. The operation conditions are: 
the rotation speed of the engine is -120.0 rpm (clockwise): the 
cold-end and hot-end temperatures are 300.0 K and 400.0 K, 
respectively. The first step of any CFD simulation is to determine 
a proper grid and time step interval to obtain grid and time-step 
independent solutions. Three different grids with cell numbers of 
80,566, 110,569, and 141,438, respectively were tested by using 
three different time-step intervals of 2.5 x 10~ 3 s, 1.667 x 10~ 3 s, 
and 1.25 x10~ 3 s, which respectively correspond to time-step 
numbers per cycle of 200, 300, and 400. Fig. 2 shows the grid with 
80,566 cells. As seen, grid density has been concentrated towards 
all solid boundaries to resolve the flow and thermal boundary lay¬ 
ers. For each simulation, about 4 cycles were needed for the solu¬ 
tion to become periodical. Judging from the values of maximum 
velocity, heat input, and work output within a cycle, grid and 
time-step independent solutions can be achieved by the grid with 
110,569 cells at 200 time steps per cycle. The three quantities 
obtained by this grid and time-step interval are less than 0.5% in 
difference with those obtained with the finest grid and smallest 
time-step interval; therefore this grid and the time-step interval 
have been used in this study. To validate the correctness of 


0=36.0° 6=72.0° 6=108.0° 0=144.0° 6=180.0° 




Fig. 5. Temperature contours across the whole engine at 10 different crank angles. 
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USTREAM code, an adiabatic cycle of the current engine has been 
performed, and the p aV e-V diagram is shown in Fig. 3. For an adia¬ 
batic cycle, the exact relation between p and V is pV 14 = C, where C 
is a constant. Fig. 3 indicates that USTREAM code returns results 
that fit this relation very well; hence, the correctness of the code 
has been validated. 

Fig. 4 shows the variations of total, expansion chamber, and 
compression chamber volumes versus crank angle over an engine 
cycle. The maximum and minimum values of the total volume 
are respectively 1.586 x 10~ 4 m 3 and 1.979 x 10~ 4 m 3 , resulting 
in a compression ratio of 1.248, and they occur at 9 = 90.0° and 
270.0°, respectively. Since the expansion chamber includes the 
power cylinders, the maximum volume of the compression cham¬ 
ber is larger than that of the expansion chamber. 

To demonstrate the general variations of temperature within 
the engine over an engine cycle, Fig. 5 shows the temperature con¬ 
tours of the entire engine at 10 different crank angles, from 
0 = 36.0° to 360.0° with an increment of 36.0°. The motions of 


displacer and power pistons pose different effects on the flow in¬ 
side the engine. The combined motions of displacer and power pis¬ 
tons cause gas to flow forwards and backwards between the 
expansion and compression chambers, while the motion of the 
power pistons expands or compresses the total volume of the en¬ 
gine, which in turn introduces cooling or heating effects on the 
working gas. The gas that flows forwards and backwards between 
expansion and compression chambers introduces very strong heat 
convection effect in both chambers and results in complicated heat 
transfer behaviors between the gas and the solid boundaries and 
highly non-uniform temperature distributions across the chambers 
at any given moment. In addition, due to the fact displacer and 
power cylinders are eccentric, the temperature distributions are 
not axisymmetric, and there exist pronounced 3D phenomena 
throughout the entire engine cycle. To this end, it is readily realized 
that the assumption of uniform temperature distribution and the 
adoption of constant convective heat transfer coefficients to ac¬ 
count for the heat transfer effect in expansion and compression 




Fig. 6. Temperature contours, velocity vectors, and local heat flux distributions in the expansion chamber at 0 = 90.0°; (a) temperature contours and velocity vectors, (b) local 
heat flux distributions at bottom plate. 
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chambers in some second-order models is apparently too crude to 
represent the real heat transfer behaviors inside the engine. There¬ 
fore, the attention is directed towards the investigation on the en¬ 
gine’s heat transfer characteristics which can only be analyzed in 
detail by using a CFD approach. 

From the expansion and compression chambers’ point of view, 
both go through one phase of receiving gas injection from the other 
chamber and the other phase of gas ejection to the other chamber 
through the regenerative channel. In the following, these phases 
are termed “injection phase” and “ejection phase”, respectively. 
In general, the injection phase of the expansion chamber is the 
ejection phase of the compression chamber, and vice versa. In 
the following, the discussion on the flow and heat transfer charac¬ 
teristics is divided into injection phase, ejection phase, and in 
regenerative channel subsections. 

4.1. Injection phase 

The injection phase in the individual chamber generally coin¬ 
cides with the expansion process of its chamber volume. In the 
expansion chamber, this is from 0 = 0.0° to 180.0°; and in the com¬ 
pression chamber, it is from 0 = 156.6° to 347.4°. As shown in Fig. 5, 
this phase in the individual chamber can be characterized by the 
appearance of a major thermal plume in that chamber. The injec¬ 
tion phase poses significant impact on the heat transfer in both 
chambers. To demonstrate this effect, Fig. 6 shows the enlarged 
views of velocity vectors and temperature contours as well as local 


heat flux distributions in the expansion chamber at 0 = 90.0°. 
Fig. 6(a) show that gas injection produces a jet of cold gas that im¬ 
pinges on the bottom heated wall of the chamber (the hot-end 
wall). This impingement brings the cold gas in direct contact with 
the hot-end wall, resulting in impingement heat transfer which is 
well known to be the most effective convective heat transfer mech¬ 
anism. Therefore, the cold gas get injected into the expansion 
chamber is rapidly heated up. This can be verified from the distri¬ 
butions of the local heat flux per unit area at 0 = 90.0° given in 
Fig. 6(b). The local heat flux q is calculated by: 


In this equation, the heat transfer from solid wall (surrounding) 
to gas (thermodynamic system) is positive, and that from gas to 
wall is negative. With referencing to the velocity vectors in 
Fig. 6(a), pronounced elevation of local heat flux can be seen in 
Fig. 6(b) at the location where the cold gas impinges on the hot 
bottom wall. The maximum local heat flux on the hot-end wall 
can be as high as 5.4 kWm~ 2 . Fig. 6(b) also shows that large local 
heat flux only occurs at the near vicinity of the impingement re¬ 
gion, and the magnitude of local heat flux quickly decreases to¬ 
wards the center of the bottom plate. The plume of cold gas 
forms a primary vortex which in turn induces some very small sec¬ 
ondary vortices at its sides. These secondary vortices promote mix¬ 
ing of cold and hot gas and further facilitate convective heat 
transfer in the expansion chamber. However, their effects are much 



Fig. 7. Temperature contours, velocity vectors, and local heat flux distributions in the compression chamber at 0 = 270.0°; (a) temperature contours and velocity vectors, (b) 
local heat flux distributions at top plate and power cylinder wall. 
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smaller compared to that of the impingement. The local heat flux 
distributions indicate that the dominant heat transfer mechanism 
in the injection phase of the expansion chamber is the impinge¬ 
ment of the cold gas on the hot-end wall. 

Fig. 7 shows the enlarged views of velocity vectors, temperature 
contours, and local heat flux distributions on the top plate of the 
compression chamber at 0 = 270.0°. Similar impingement phenom¬ 
enon can be observed during the injection phase of the compres¬ 
sion chamber, but the impingement results in rapid cooling (local 
heat flux is negative) of the hot gas that enters the compression 
chamber. Here the maximum absolute local heat flux at the top 
plate during the injection phase is 3.0 kWm~ 2 . This is much smaller 
than that in the expansion chamber during its injection phase. 
However, due to the compression effect posed by the movement 
of power pistons, there is large negative local heat flux occurring 
at the ejection phase, which will be discussed in the next sub-sec¬ 
tion. Here, the impingement heat transfer remains the dominant 
heat-transfer mechanism in the compression chamber during its 
injection phase. 

4.2. Ejection phase 

The ejection of gas in expansion or compression chambers is 
mainly caused by the reduction of volume in the individual cham¬ 
ber. By referring to Fig. 4, this phase is from 6 = 180.0° to 360.0° in 
the expansion chamber and from 0 = 347.5° to 156.5° in the com¬ 
pression chamber. As mention earlier, the reduction of chamber 
volume is due to the motions of the displacer and power pistons. 


When displacer moves towards the hot end, the volume of expan¬ 
sion chamber is reduced. On the other hand, when displacer and 
power pistons move toward each other, the volume of compression 
chamber is reduced. In this manner, the displacer and power pis¬ 
tons reduce the volume in the individual chamber by decreasing 
the distance between the upper and lower boundaries of that 
chamber, and this is equivalent to push the gas against the solid 
walls. In the expansion chamber, the lower wall is the hot end, 
while in the compression chamber; the cold end includes the 
upper wall and the walls of power cylinders. The pushing of gas 
by the displacer and power pistons creates an effect similar to 
impingement that largely compresses the isothermal lines towards 
the hot or cold ends and increases the local temperature gradient. 
As a result, it also promotes heat transfer in a way similar to flow 
impingement does. It can be seen in Fig. 4 that in the compression 
chamber, the ejection phase mostly coincides with the compres¬ 
sion process of the entire engine volume, and the negative heat 
transfer is further boosted by the compression effect which heats 
up the gas in the compression chamber and creates even larger 
temperature difference between gas and cold-end wall. The com¬ 
bined heat transfer effects can be demonstrated by velocity vec¬ 
tors, temperature contours, and local heat flux distributions in 
the compression chamber at 0 = 54.0° shown in Fig. 8. The temper¬ 
ature contours in Fig. 8(a) show that the core temperature in the 
compression chamber reaches 330.0 K due the compression pro¬ 
cess. Consequently, the negative local heat flux on top plate and 
wall of power cylinder is largely boosted as can be seen in 
Fig. 8(b). The absolute maximum heat flux reaches 6.8 kWm 2 , 




Fig. 8. Temperature contours, velocity vectors, and local heat flux distributions in the compression chamber at 0 = 54.0°; (a) temperature contours and velocity vectors, (b) 
local heat flux distributions at top plate and power cylinder wall. 





















































W.-L. Chen et al. / International Journal of Heat and Mass Transfer 75 (2014) 145-155 


153 


which is more than twice the maximum value during the injection 
phase. In the expansion chamber, on the other hand, the ejection 
phase does not coincide so well with the expansion process of 
the entire engine volume; therefore, the promotion on positive 
heat transfer by the mechanism of enlarging temperature differ¬ 
ence between gas and hot-end wall via gas expansion is not very 
significant. On the contrary, Fig. 4 shows that from 0 = 269.0° to 
360.0°, the ejection phase of the expansion chamber coincides with 
the compression of total engine volume, thus the gas inside the 
expansion chamber is heated not only by the hot-end wall but also 
by the compression effect. As a result, the local maximum temper¬ 
ature in the expansion chamber reaches as high as 417.0 K, which 
is higher than the hot-end temperature. The elevated gas temper¬ 
ature even gives rise to small negative heat flux on bottom 
plate. Due to the above factors, the local heat flux on the bottom 
plate is moderate during the ejection phase of the expansion 
chamber. 

4.3. In regenerative channel 

The major function of a regenerative channel is to pre-heat the 
gas during the injection phase of the expansion chamber and to 
pre-cool the gas during the injection phase of the compression 
chamber. From the local heat flux distributions on the wall of dis¬ 
placer cylinder shown in Fig. 9, it can be seen that the regenerative 
channel performs this task relatively well as the heat flux during 
the injection phase of expansion chamber on regenerative wall 
shown in Fig. 9(a) is mostly positive (heating gas), and that during 
the injection phase of compression chamber shown in Fig. 9(b) is 
mostly negative (cooling gas). However, there are exceptions, in 
terms of the injection phase of the expansion chamber, Fig. 9(a) 
shows that negative heat flux is present at the upper section of 
the displacer cylinder wall. This is because gas at a higher temper¬ 
ature than the wall temperature entering the regenerative channel. 
The reason for the elevated gas temperature is due to the heating 
from the compression of the entire engine volume. In terms of 
the ejection phase of the compression chamber, Fig. 9(b) shows 
that at the lower section of displacer cylinder wall (near entrance 
of the regenerative channel), heat flux is positive. This happens 
when the gas with a temperature lower than the wall temperature 
in the expansion chamber entering the regenerative channel. The 
lower gas temperature in the expansion chamber is due to the 
cooling of gas by the expansion of the entire engine volume. Due 
to the eccentric location of the power cylinders, Fig. 9 again shows 
that the local heat flux distribution is not axisymmetric, especially 
near both entrances of the regenerative channel. 

4.4. Average and integrated quantities 

The p aV e-V diagrams of the expansion and compression cham¬ 
bers are illustrated in Fig. 10. Here p ave is the averaged pressure 
in the individual chamber. It can be seen that during one engine cy¬ 
cle, the expansion chamber yields positive work, whereas the com¬ 
pression chamber produces negative work. The absolute value of 
the positive work is larger than that of the negative work, giving 
rise to positive net work output to the engine’s shaft. The varia¬ 
tions of average temperature in the expansion and compression 
chambers are shown in Fig. 11(a), while those of heat flux are 
shown in Fig. 11(b). Here heat transfer rate Q_ is calculated by 
integrating the local heat flux q over the individual chamber’s wall 
surface. Fig. 11(a) indicates that average temperature fluctuation in 
the expansion chamber is within the range of 28 K, and due to the 
compression effect mentioned earlier, the maximum average tem¬ 
perature is 401.0 K which is slightly higher than the hot-end tem¬ 
perature. In the compression chamber, the average temperature 
fluctuation is within the rage of 20 K. Both average temperature 


(a) 


e= 90.0° 



Fig. 9. Local heat flux distributions on the regenerative wall at two different crank 
angles; (a) 6 = 90.0°, (b) 9 = 252.0°. 


variations show only one peak and one valley, corresponding to 
one injection phase and one ejection phase within an engine cycle. 
In Fig. 11(b), it can be seen that heat transfer rate in the expansion 
chamber in the injection phase first increases steadily, reaches the 
maximum, and then decreases monotonously. In the ejection 
phase, however, the general tendency is decreasing but it briefly 
rebounds twice during the descent. This is due to the fact that 
the second half of the ejection phase of the expansion chamber 
generally coincides with the expansion of the entire engine vol¬ 
ume, and the expansion effect cools down the gas temperature, 
thus promoting positive heat transfer. The amount of heat transfer 
rate (in terms of absolute value) in the compression chamber also 
shows monotonous increasing then decreasing in its injection 
phase, in its ejection phase, the general tendency is decreasing 
but there it also briefly rebounds once due to a period of coinci¬ 
dence of the ejection phase with the compression of the entire en¬ 
gine volume as described earlier in the Section 4.3. 
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Due to the expansion and compression of engine volume, heat 
input and output can happen at the hot end, cold end, and on the 
wall of regenerative channel. Therefore, the integration of heat in¬ 
put cannot be confined just along the wall of the expansion cham¬ 
ber. Neither can the integration of heat output be confined along 
the wall of compression chamber. The rates of heat input and out¬ 
put are calculated by integration along all solid boundaries: 

to /' t=t » +t p r 

Q.in = T~ / qdAdt, if g > 0, (18) 

Jt=t 0 Jwall 


m ft-to+tp n 

Q out = — / / qdAdt , if g < 0. (19) 

In terms of the evaluation of output power, the focus is on the 
surface of the power piston because this is the exact place engine 
work is delivered. The power is calculated by: 

where the negative sign on the right hand side of equation is to ac¬ 
count for the negative to in this study, and the circular integration is 
from minimum z p to maximum z p and back to minimum z p . The cal¬ 
culated rates of heat input, heat output, and power are 13.28 W, 
12.63 W, and 0.64 W, respectively, resulting in a thermal efficiency 
q of 4.82%. Here, q is calculated by dividing the engine power by the 
rate of heat input. Such a low efficiency is mainly due to the small 
temperature difference between the hot and cold ends and the low 
engine speed. Since the focus of the study is to understand the heat 
transfer characteristics of the Stirling engine cycle, to improve the 
engine’s efficiency will be dealt with in some future study. 

To help readers understand the dynamic processes inside the 
current Stirling engine, two videos have been made. The video 
named "gamma_aH”, plays the transient variations of temperature 
contours and velocity vectors in expansion and compression cham¬ 
bers within two engine cycles. The other video named “flux_all” 
plays the variations of local flux on bottom plate, top plate, power 
cylinder wall, and regenerative wall, within two engine cycles. 




eft 


Fig. 11. Variations of average temperatures and heat transfer rates in the expansion 
and compression chambers; (a) variations of average temperatures, (b) variations of 
heat transfer rates. 


5. Conclusion remarks 

An in-house CFD code has been developed and applied to a 
twin-power-piston y-type Stirling engine to study the heat transfer 
characteristics in the thermodynamic cycle of the Stirling engine. 
This is the first CFD simulation on this engine, and the heat transfer 
characteristics have been investigated in great detail. The results 
include temperature contours, velocity vectors, and distributions 
of local heat flux at several important stages in an engine cycle 
as well as average temperature variations, integrated rates of heat 
input, heat output, and engine power. Some conclusion remarks 
from this study can be drawn as follows: 

1. The complicated heat transfer behaviors in the engine result in 
highly non-uniform temperature distributions in all sections of 
the engine at most of the time in an engine cycle. The 
assumption of uniform chamber temperatures and the use of 
constant heat transfer coefficients to account for heat transfer 
adopted by some zero- or one-dimensional models are too 
crude to reflect the reality. 
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2. In the injection phase, impingement heat transfer is the major 
heat transfer mechanism which allows gas to be heated (in 
expansion chamber) or be cooled (in the compression chamber) 
very effectively. However, this effect tends to concentrated at 
some localized region, resulting in non-uniform heating or cool¬ 
ing of gas. 

3. In the ejection phase, the motions of displacer and power piston 
also create a sort of heat transfer effect similar to impingement, 
but this effect acts on a much larger region and also promotes 
heat transfer effectively especially in the compression chamber 
where the negative heat transfer is further boosted by the rising 
of gas temperature due to compression of the entire engine 
volume. 

4. The regenerative channel functions as it was designed for at 
most of the time within an engine cycle. It pre-heats or 
pre-cools gas before gas got injected into the expansion or 
compression chamber. However, exceptions can occur near both 
regenerative channel entrance regions at some brief moments. 
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